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Abstract 

For Major ana- Wilson lattice fermions in two dimensions we derive a 
dimer representation. This is equivalent to Gattringer's loop representa- 
tion, but is made exact here on the torus. A subsequent dual mapping leads 
to yet another representation in which a highly efficient Swendsen-Wang 
type cluster algorithm is constructed. It includes the possibility of fluctuat- 
ing boundary conditions. It also allows for improved estimators and makes 
interesting new observables accessible to Monte Carlo. The algorithm is 
compatible with the Gross-Neveu as well as an additional Z(2) gauge inter- 
action. In this article numerical demonstrations are reported for critical free 
fermions. 
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1 Introduction 



Occasionally in talks or papers about dynamical fermions it is mentioned — more 
or less as a joke — that the computer has no data-type Grassmann and one hence 
can simulate fermions only via the nonlocal effective theory after integrating them 
into the determinant. Of course, this is plagued by the well-known inefficiencies. 
In this article, based on Gattringer's loop representation [1], we show that in two 
space-time dimensions one actually can get pretty close to 'simulating Grassmann 
numbersjl]. We here expand on py and its recent numerical implementation [2] in 
several ways. First the loop representation is re-derived starting from Majorana 
fermions in what we think is a particularly natural way. The new connection in- 
cludes definite boundary conditions on the torus and does not only work in the 
thermodynamic limit as before. In particular, we can then also approach the fi- 
nite volume continuum limit. Furthermore we propose a cluster algorithm that is 
(practically) free of critical slowing down and allows for improved estimators. In 
this formulation we can simulate fluctuating boundary conditions which is neces- 
sary to allow for fixed (anti)periodic boundary conditions in the original fermion 
system. It also makes ratios of partition functions accessible as observables in 
Monte Carlo simulations. They constitute interesting quantities in the continuum 
limit. 

The original Gross-Neveu model of self-coupled fermions in two dimensions [S] 
is most naturally written in terms of N species of Majorana fermions. In the lattice 
discretization with Wilson fermions the euclidean action is given by [1] 

^ = E { l^'^^M + ^ - '2^9*8)^ - ^jiCco'] ■ (1.1) 

The Grassmann- valued field ^ = ^01(2;) has a spin index a = 1,2 and a flavor index 
i = 1, . . . N that we leave implicit. We denote by d,d*, d the forward, backward 
and symmetric nearest neighbor differences on our cubic T x L lattice. The charge 
conjugation matrix C obeys 

C7,C-^ = -7j = -7;, C = -C^. (1.2) 
For even each pair of Majorana fermions may be considered as one Dirac fermion 

with its independent In the Majorana form the full global symmetry group 

0{N) is manifest beside (without Wilson term) the discrete 75 symmetry whose 

^Cum grano salis. 
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spontaneous breaking was studied in [3] in the N ^ oo limit. The model is 
renormalizable in the strict sense: there is no other O(A^) invariant scalar 4-fermion 
interaction term. For iV = 2 we have the Thirring model [5J, [6J, the cases N ^ 3 
are expected to be asymptotically free. In the remainder of this paper we set the 
Wilson parameter to r = 1 and work in lattice units a = 1. The discrete chiral 
symmetry ^ 75^ of the massless continuum theory is broken by the Wilson term 
and is only expected to be recovered in the continuum limit at the critical mass 
m = rric- On the torus we consider four conceivable combinations of periodic or 
antiperiodic boundary conditions in the two directions. Periodicity angles different 
from 0, TT — as sometimes used for Dirac fermions — would not lead to a periodic 
action density for Majorana fermions. We label the possible boundary conditions 
by a bit-vector 

e^, eo,eie{0,l} (1.4) 

where stands for periodic and 1 for antiperiodic boundary conditions in the 
corresponding direction. 

Often the interaction term is factorized by the introduction of an auxiliary 
bosonic field. For us it will be more convenient to think of m — m{x) as an x- 
dependent mass for a while. If Z|[m] is the partition function of one free Majorana 
fermion with boundary condition in this background field, then the partition 
function of the interacting theory is written as 



2 dm{x 

Integrating the fermions in the remaining Gaussian problem yields the Pfaffian 



Z|[m] = Pf 



C{^A + m-^d*d) 



;i.6) 



In appendix |A] one can find a reminder of the definition of Pfaffians. For even 
we may replace the Pfaffians by the N/2 power of the determinant. In this form 
and with a factorizing field, the model can be simulated by standard methods 
like HMC as carried out in [1]. For larger g this compute- intensive task became 
rather difficult due to singularities developing in the operator under the Pfaffian. 
Of course, the model itself remains completely well-defined on any finite lattice. 
Fermionic and compact bosonic variables are safe in this respect. 

As a first step we now introduce the loop representation [1] of Z|[m]. It may 
in fact also be looked upon as a dimer ensemble similar to those derived in [7] for 
strong coupling QCD. 
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2 External field fermion partition function 

As a building block for the Gross-Neveu models we consider for a single Majorana 
field the external field action 

s = ^Yl vix)e{x)mx) - J2 ii)e{x)cp{m{x + a)- (2.1) 

X x,ix 

We assume a lattice with T sites in the time direction {fi = 0) and L sites in 
the space direction (/x = 1). The variables (p{x) = 2 + m{x) and r(x,/i) are 
external commuting fields. The link field r is introduced for completeness. It 
will be dropped at some point. The lattice derivatives in (11.11) combine to Wilson 
projectors that we define for arbitrary lattice unit vectors n = ±fi 

Pin) = i(l - n,^,), {CP{n)Y = -CP{-n). (2.2) 

The last identity implies 

e{x)cpm{x +jj) = e{x+ fi)cpi-mx)- (2.3) 

While the field ^ has torus periodicity e, the external fields (p, r are continued 
periodically to obtain a periodic action density. 
Defining the 'covariant' projecting hop operator 

(H^Oix) = rix,fi)CP{mx + fi) (2.4) 

we may, with the help of (12. 3p . also write the action in the manifestly antisymmetric 
short- hand form 

^ = ^ E V'^^^^ - ^ E ^^(^M - Hj)^ (2.5) 

X X,fl 

where Hj is transposed with respect to both spin and space indices yielding 

(HjO (x) = -rix - A, /i)CP(-/i)e(a; - /i). (2.6) 
Now we are ready for the partition function 

ZI[ip,T]= j Die-'. (2.7) 

The measure is 

Di = \[{di^di2){x) (2.8) 

X 
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and yields 

Z^[cp,T] = Fi\cpC-y^{H,-H;,) I (2.9) 



which is a nonlocal expression in the external fields as in the case of the usual 
Dirac fermion determinant. In appendix El we evaluate the Pfaffian for r = 1 and 
constant for all four boundary conditions. 

3 Equivalent statistical systems 

3.1 Dimer representation 

The Grassmannian Boltzmann factor may be expanded, 

[ D^i[{i+^^2^,}i[{i+T{x,f^)e{x)cpim^+f^))- (3.1) 

X X,fJ, 

All fields in the curly bracket are at x and this factor is best considered as part of 
the measure. We have here chosen C such that \^~^C^ = ^1^2 which just amounts 
to a phase convention. Note that the square of the hop-term vanishes due to 
the one-dimensional projectors. There is only one linear combination of the two 
Grassmann numbers contributing from each site which squares to zero. We now 
introduce one-bit-valued dimer or bond variables [7] on each link k{x,fi) = 0, 1, 
whose values are used to organize the expansion as 



{k(x,^i)} X x,fi 



(3.2) 

As in [7J the goal now is to integrate out the fermions to yield a Boltzmann weight 
p[k] for each dimer configuration. By asking how the Grassmann integrations can 
be saturated site by site it is clear that a non-zero weight only arises if at each 
site there are either two dimers adjacent from different links or none at all. In 
the latter case the integration is saturated by the measure term and a factor <f{x) 
appears for this site. We also call these contributions monomers. Due to the above 
constraint the dimers have to form closed non-intersecting and non-backtracking 
loops. We choose randomly a starting point and an orientation on such a loop such 
that along the loop one visits the sites (xi, X2, ■ ■ ■ , xi). Consecutive sites differ by 
lattice unit vectors, Xj+i = Xi + hi, including xi = xi + hi in the last step. For 
such a loop the product of bilinears 

{e{x,)CP{h,)i{x,)) {C{x2)CP{h2)ax,)) ■ ■ ■ {C{xi)CP{hi)ax,)) (3.3) 
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has to be considered together with the integrations on the / sites involved. Note 
that here (12. 3p is relevant on the links that are transversed in the negative direction. 
The trivial key formula is, for a single site, 

j d^id^2^e = C-\ (3.4) 

Now the expression (13. 3p integrates to 

X = -tT[Pihi)P{h2)---Pini)]. (3.5) 

Here appears the very important minus sign for a closed fermion loop, well-known 
for instance from Feynman diagrams. It is not difficult to see that this result is 
independent of the starting point and orientation chosen. 

The evaluation of the spin factor follows [H]. Let us introduce eigenspinors of 
the projectors 

P{ni) = \ni) (nil, {hi\hi) = 1, (3.6) 

X = -{ni\h2) {n2\n3) ■ ■ ■ {ni\hi) (3.7) 

A spinor is rotated by an angle 6 by the unitary spin matrix R{6) = exp[(0/2)7o7i]. 
This allows us to write 

\hi) = R{A9i)\hi+i) A9,e{0,±7i/2}, z = l,2,...,l (3.8) 

with ri/+i = ni. Using 

{hj\R{A9i)\nj) = cos{A9i/2) (3.9) 

we get to 

X = -(n/|ni) JJcos(A^,/2). (3.10) 

1=1 

The rotation accumulated in steps (13.81) is 

\hi) = R{Q - A9i)\hi) with e = ^A^,. (3.11) 

1=1 

For closed paths we have G = 27rz/ and R{Q — A9i) = cos{tti')R{—A9i). For the 
nonzero lattice angles, cos(A^j/2) = 1/V2 . Altogether the final result is 

X = (_)-+i2-^=/^ (3.12) 

where is the number of ±7r/2 angles ('corners') occurring along the loop and 
= 0, ±1 is the number of complete rotations the loop makes. The extra minus 
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sign for |z/| = 1 is the one associated with fermions under 27r-rotations. If we 
include a non-trivial r field then the product of its Wilson loops for all dimer 
loops appear in addition in the weight. After this remark we set r = 1 until 
further notice. 

Although we are on a lattice here, we can define homotopy classes of loops. 
Two loops are homotopic to each other if they can be transformed into each other 
by a sequence of steps, where dimers are only changed around a single plaquette. 
We see from the above that X is positive for all configurations containing only 
loops that are homotopic to the trivial loop, just a point. The two minus signs 
characteristic of fermions compensate each other in this class! 

Loops can wind however around the torus in either direction as noted in [2]. 
A pair of loops winding around the same direction is still in the trivial homotopy 
class. An odd number of windings leads however to a new class. This may also 
happen in both directions at the same time and hence there are the four classes 
£oo; -^10; '^01; -^ii- flgurc [T] wc show a representative of each class. They are 
equilibrium configurations of free fermions at m = and T = L = 10. The meaning 
of the +/— signs in the plots will become clear later. Only configurations from £oo 
have a positive weight while in the other cases there is an odd number of closed 
loops with zero total rotation angle each of which contributes a factor —1. By 
introducing antiperiodic boundary conditions in some direction the loops closing 
around that direction receive yet another sign without changing the topologically 
trivial ones. 

With the local weight 




(3.13) 



X 



with 




ip{x) if monomer at x 

I if 2 dimers in the same direction at x 

I I \/2 if 2 dimers in different directions at x 
else 



(3.14) 



we now define the positive dimer partition functions 





J2 P{k], etc. (3.15) 



{k{x,fj.)}eCoo 



From what was said above the connections 





(3.16) 
(3.17) 




(3.19) 
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Figure 1: Examples of dimer configurations (solid lines) in the topological classes 
£oO) -^10; -^oi; -^11 (from Upper left to lower right). The time direction (/i = 0) 
points to the right. The signs refer to spins on the dual lattice. 



arise, which can be inverted. If we want to realize the boundary conditions = 
(1, 0) of [4J (or actually any other definite choice for the fermions) we have to sum 
over all dimer classes including negative weight contributions 



> [^] = Iv] + {^\ - Zf {^\ + Zl' {^\ . (3.20) 

All these relations between partition functions can be turned into relations 
between expectation values of the scalar fermion density and monomer densities 
by differentiating with respect to ^{x). One example based on (13.201) is 



(K(x)) 



00 



(K(x)) 



10 



{K{x)) 



11 



+ zi^ 



Zf + 



(3.21) 
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The observable K{x) is one if there is a monomer at x and zero otherwise. With 
the help of r(a:, n) as a source one could establish further relations. 

For free fermions {ip = 2 + m) in the thermodynamic limit at fixed m > the 
various Z| differ only by exponentially small amounts and dominates among 
the dimer ensembles. Taking the finite volume continuum limit (L ^ oo with Lm 
fixed, see appendix [X]) and in particular for m = this is not so. In the latter 
case we have an exact zero fermion mode for = (0, 0) and = holds. This 
implies 

= + + Zl^ for m = 0. (3.22) 
3.2 Spin representation 

In this subsection we transform the dimer system to yet another representation by 
Ising spins. This will allow us to design a global cluster algorithm. A clue that this 
may be possible is given by the idea that a natural way to manage and modify the 
closed loops in the dimer formulation is to consider them as boundaries of domains 
of up-spins surrounded by down-spins (Peierls contours). 

The spins that we introduce live on the lattice dual to the one carrying the 
fermionic and the dimer variables. Its sites labelled by underlined x_ are dual to 
plaquettes of the original lattice and are imagined to be located at their centers. 
Analogously, the sites of the original lattice are dual to plaquettes in the new one. 
Links (x, /i) and (x, /i) are dual to each other if they cross, see figure O The idea 
is now to put an Ising field s{x) on the dual lattice and to identify configurations 

k(xu)-[^ s(x)s(x + ^) = -l 

s(x)s(x + ^) = +l • ^'^■^'^> 

In other words, dimers are located where nearest neighbor spins on the dual lattice 
are antiparallel. In a first stage we restrict ourselves to the class £oo of dimer 
configurations. 

We first prove that for each admissible dimer configuration there are exactly 
two spin fields obeying fl3.23p that differ by a global spin-flip. In a first step we 
define a Z(2) lattice gauge fielci on the dual lattice in terms of the dimers on the 
original lattice 

;[ (3.24) 

Because of the constraints on k this gauge field is unity when multiplied around 
any plaquette (on the dual lattice). As we restrict ourselves to £00 here, also loops 
around the torus are unity for this gauge field. Thus cr it a pure gauge on the 

■^This field has nothing to do with the earlier r. 
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Figure 2: Illustration of the labelling of the original and the dual lattice. 



torus with a periodic gauge function. The spin-field is this gauge function that we 
construct now. We choose a site y (the origin, for example) and set s{y) = +1. 
Now this value is parallel-transported with a to all other sites, for instance along 
a maximal tree rooted at y. Due to the absence of curvature in a, the result is 
path- independent, consistent and unique. Starting from s{y) = —1 we obtain the 
other configuration associated with {k{x,fi)}. The signs in figure [T] are just these 
spins. 

While we now have exactly two spin fields for each admissible dimer configura- 
tion in Cqq, not all conceivable spin fields are reached in this way. Obviously, spin 
configurations, that on any plaquette look like 

do not occur in the image, where we here just gave simple labels to the four spins. 
They would correspond to crossing loops not allowed by the original Grassmann 
variables. If this is however excluded on all plaquettes, then we can reconstruct 
admissible {k{x,fi)} configurations from {s{x)}. The total Boltzmann factor in 
the spin representation is now a big product with one factor for each plaquette. 
These weights w are given in table [1] which lists only 8 of the 16 configurations 
and is completed by using w{si, S2, S3, S4) = w{—Si, —S2, —S3, — S4). The monomer 
strength (f{x) is taken here at the site x of the original lattice sitting at the center 
of the dual plaquette considered. 

The derivation of this representation resembles the construction of the dual 
formulation for generalized Ising models p]. We summarize it (giving the ordi- 
nary self-dual two-dimensional Ising model as an example in brackets): One first 
introduces new variables living on the bonds making up the Hamiltonian of the 
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Table 1: Plaquette weights depending on the spin configuration. 



original theory (a link field). Then the original spins are summed over producing 
a constraint in the new variables (vanishing plaquettes of the links interpreted on 
the dual lattice). This constraint is then solved on the dual lattice (links given 
as a pure gauge by a site field). The extension of the concept here is that we 
change from Grassmann elements to bosonic variables, have an additional con- 
straint fl3.25p to fulfill, and that there can be minus signs. One could talk of 
Fermi-Bose- or super-duality. 

The plaquette weight can also be written in terms of pairwise nearest neighbor 
bond-interactions of the form [writing 5ij = 6s^sj = |(1 + SiSj)] 

yj{Si,S2,S3,Si) = {p[5l2 + §23 + + hi] + q[Sl2S34 + Su523\ + 

r[6i26u + ^21^23 + ^32^34 + ^41^43]} , (3.26) 

with the x-dependence in p, q, r suppressed. To match table [Tj the following equa- 
tions have to hold 

ip{x) = Ap + 2q + 4r 
1 

—= = 2p + r 
1 = 2p + q. 

The solution of this system is 

1 TTl \ T 1 

r = -((/?- 2) = —, p= — q = l ^ + r. (3.27) 

4^^ ' 4' 2^2 2' ^ 72 

We remark here that all coefficients are positive for ^ m ^ 2^2 ^ 2.83. In the 
free case, this clearly covers the relevant range of bare masses. 

If we now include dimer configurations in the other sectors by using fl3.23p . 
the only difference is that the resulting s{x) is antiperiodic in the direction or- 
thogonal to those where dimer loops run around the torus. Thus Ciq corresponds 
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to spin-fields antiperiodic in space, £01 to those in time and £11 to spins with 
both directions antiperiodic. The mechanism here is that antiperiodic boundary 
conditions of the spins force an interface into their configurations which leads to 
the nontrivial dimer loop topology. Introducing further partition functions for 
the spin ensembles we relate 

^fc — 2 s ' I' ~ 2 " ' ^ ~ 2 ' ~ 2 ' l^-^oJ 

The factor 1/2 cancels the global spin-flip symmetry. Again, by differentiation, we 
may relate expectation values, where the presence (absence) of a monomer yielding 
K{x) = 1(0) in the spin language translates into maximally 'polarized' plaquettes 
where all four spins are parallel, 

K{x) =6^,\Mix)\, M{x)= ^(^)- (3-29) 

X around x 

4 Simulation algorithms 
4.1 Local algorithms 

A local algorithm to simulate any one of the above dimer ensembles with all weights 
taken positive was recently described and tested by Gattringer et al. [2j. The 
simplest case to consider is a free Majorana fermion of mass m with (p{x) = 2 + m. 
In the updates one actually performs only changes that are local in the homotopy- 
sense by proposing dimer-fiips k{x, /i) 1 — k{x, fi) around plaquettes. The move 
is accepted with the Metropolis probability corresponding to the ratio of p in fl3.13p 
for the new and the old configuration, which is a locally computable quantity. Of 
course, it vanishes, if the new configuration would violate one of the constraints. 
This update stays in the homotopy class fixed by the starting configuration (see 
figure [T]) and is ergodic within it. Thus the various ensembles can be simulated 
which correspond to combinations of periodic and antiperiodic Pfaffians. In [2] 
it was demonstrated that already such simulations are vastly more efficient than 
HMC type simulations. 

The entirely equivalent update in the Ising form consists of local spin-flips. The 
Metropolis decision in this case depends on the eight nearest and next-to-nearest 
(diagonal) neighbors that share plaquettes with the spin in focus. The numerical 
efficiency in both forms is very similar. 
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4.2 Cluster algorithm 



The plaquette interaction in f l3.26p is now written as a superposition of 10 different 
terms, schematically given by 

10 

w{si,S2,S3,Si) = ^PiAi{si, 82,33,84) (4.1) 
i=l 

with Pi G {p,q,r} and A(si, S2, S3, S4) G {0,1}. In complete analogy to we 
introduce now ten-valued bond variables b{x) with each value corresponding to 
one of these terms to obtain 

Zs= ^ ]^-P6(x)Afe(^)(si,S2, 53,34). (4.2) 

{b{x), s{x)} plaq 

To avoid a too clumsy notation here again 51,52,53,54 are spins around each of 
the plaquettes. The celebrated trick of [ID] consists of Monte Carlo sampling both 
the b and the 5 variables. As the first part of an update cycle one chooses new 
b at fixed 5 by a local heatbath procedure. Of course, in general, the choice is 
between less than 10 possibilities, as some of the Ab(j,.)(5i, 52, 53, 54) vanish. Then, 
for given bonds b, any of the 2'^^ spin configurations has either weight zero or a 
constant nonzero weight depending on the constraints given by the product of 
all factors Ab(x)- Just as in the standard Ising model case we now construct the 
percolation clusters defined by the active bonds in all {Ai,[x)}- By flipping the 
spins in each cluster as a whole with probability 1/2, we sample one of the allowed 
and equally weighted spin configurations. The overall procedure amounts to a 
global independent sampling of spins (at fixed bonds) and will be numerically 
demonstrated to almost eliminate critical slowing in section [5l 

At the stage of selecting a new spin field one may also construct improved clus- 
ter estimators. This is achieved, if, for some observable, one is able to analytically 
average over all conceivable spin-assignments of which only one is taken as the 
next configuration. 

While the above procedure is analogous to the well-known Swendsen-Wang al- 
gorithm [in] we could also study a single cluster variant [H]: one spin is chosen 
at random, and then only the one cluster connected to it is constructed by inves- 
tigating the plaquette terms touched in the growth process until it stops. Then 
spins on this cluster are always flipped. This may well be even more efficient as 
large clusters are preferred. 

We end this subsection with a remark on the global spin flip symmetry. At 
first one may think that it is a (slightly) annoying redundancy in the new repre- 
sentation. However, it is in fact essential to be able to grow clusters whose energy 
(action) is associated with the surface (in our case the loops) and not with the 
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bulk. Loosely speaking, as one grows a (single) cluster, there is always an energetic 
'way out' by flipping the whole lattice. Of course, if this is all that happens, that 
algorithm will not be efficient. The auxiliary percolation problem allows to find 
nontrivial clusters. 

4.3 Fluctuating boundary conditions 

In (13.201) we saw that it is desirable to also be able to simulate enlarged ensembles 
where one sums beside configurations of spins also over several possible boundary 
conditions. If in conventional simulations one proposes a change of the boundary 
conditions at fixed spins, one generates energy (action) proportional to T or L 
and the proposal will practically always be rejected. It was noticed however in 
|12j that with cluster algorithms the situation can be different. In the step where 
we pick new spins at fixed bonds the search among possible equally weighted new 
configurations can be enlarged to also include changed boundary conditions. If we 
label them by £^ again (for the spins now) the four possible e become a dynamical 
variable. In p!2] these changes were introduced in a single-cluster /Metropolis spirit, 
which would also be possible — in fact less involved — here. In view of the future 
construction of improved estimators we stick however to the many-cluster view 
and design now a correspondingly generalized cluster algorithm. It consists of the 
following steps: 

• We throw bonds on the links as discussed before. 

• We determine by some percolation algorithm (e.g. tree search) the indepen- 
dent spin clusters connected by bonds hut ignore two layers of links such that 
the torus is cut open. We take {{x,Q)\xq=t-i} and {(x, 1)|xi=l-i}- We call 
these clusters preclusters, their connectivity is determined in the 'interior'. 
Each of them carries a unique cluster label as a result. 

• Now the remaining links are examined as far as they have been activated. 
We call these bonds clamps. They have the effect of sewing up (some of) 
the preclusters. This is done by the pointer technique described in [I^. We 
may visualize the process as a graph with the preclusters as blobs, some of 
which get connected by lines. 

• In this process a la [13] one can detect when closed loops in the graph are 
formed. We set one of four types of flags whenever a loop is closed. They 
distinguish whether an odd or an even number of temporal or spatial clamps 
are met around the loop. We end with flags /oo, /lo, /oi; /ii each being zero 
or one. 
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• Boundary conditions e are only compatible with this loop structure if for all 
occurring loops with flag up {fap = 1) the condition aeo + [3ei = 0mod2 
holds. Of course, e = (0, 0) is always allowed. 

• Among the compatible e (1,2 or 4 values) one is chosen with equal probability. 

• Now flips for all preclusters are determined and executed. Connected com- 
ponents of the graph flip together, but preclusters within these components 
can flip relative to each other if the boundary condition has changed. The 
above construction guarantees that the orientations thus propagated do not 
depend on the path that is taken on the graph. 

Although rather short and compact in the end this is not a trivial code to write. 
It is helpful to organize it under a geometric point of view focussing on parallel 
transport between preclusters with a Z(2) group, where the boundary conditions 
are gauge variables on the clamps. Of course, as we can solve the free fermion 
ensembles exactly (appendix Ej) and easily get high accuracy with cluster simula- 
tions, many significant checks by short simulations (taking seconds) were available. 
A very good monitor for debugging at every stage is to set traps for the occurrence 
of illegal plaquettes fl3.25p . An alternative strategy to move boundary conditions 
would be to turn the loop structure of the above graph into a system of linear 
equations in the Galois field of two elements (addition isomorphic to logical xor). 
Following ^Ij and [I5] this can be solved by Gauss elimination. 

The above scheme contains some nested pointer operations. One could be 
worried in principle whether the execution time grows more than linearly with the 
lattice volume. In practice there was found to be absolutely no problem of this 
kind. This is in fact the same for cluster simulations of the standard Ising model 
using the algorithm of [13]. In our case the problem is even less severe as the 
number of clamps is smaller than T + L, not proportional to the volume. 

4.4 Negative mass 

For free fermions one could content oneself with the parameter range m ^ = rric. 
On the other hand all results in appendix |A] can be taken at arbitrary m. After all, 
the partition function on the finite lattice is just a polynomial. When we later come 
back to the interacting theory, it will also turn out, that negative mass fermions 
will be required because < due to renormalization. The local algorithm [2] 
works for negative m, a sign problem only arises if = 2 + m changes sign. The 
bond probabilities (13.271) however restrict the cluster algorithm so far to m ^ 0. 
Luckily, there is an alternative decomposition of the plaquette interaction into 
bonds that comes to rescue when m{x) is negative. 
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The r-term in (13.261) is replaced by 

-[SuSlA + 521^23 + ^32534 + 541^43 + 512^14 + 521^23 + 532^34 + 541^43], 

where we introduced also antibonds 6ij = 1 — 6ij. The two other terms are un- 
changed but the coefficients are now p, q. Using 6126^ = ^12 — ^12^14 etc. one sees 
that the new weight coincides with the old one if we identify r = —r and p = p + r. 
The matching equations are now 

ip{x) = Ap + 2q 
1 

2p + r 



%/2 

1 = 2p + g + 2r, 



solved by 



_ 1 / —TTi If 1 _ 

-=4(2-^) = ^, P=^--2^ ' = '-7^-'- 

Now all weights are positive for —m ^ 2(2 — -\/2) ~ 1.17. 

The presence of antibonds is compatible with the cluster search including fluc- 
tuating boundary conditions. With an m{x) that changes sign over one lattice, 
one actually decomposes some plaquettes with (I3.27P and others with (14. 3p . One 
must however not make the mistake to think of the preclusters as ferromagnetic 
(Weiss) domains, they contain in general both up and down spins. This is why 
we took care to talk about distributing flips to clusters rather than assigning new 
spin orientations to them as whole. 



5 Numerical applications 

We now report on numerical experiments. In this ffist publication on the method 
we stick to free fermions and the observable K{x). It corresponds to the scalar 
fermion density and is mainly used to diagnose the algorithm. Hence, with the 
results of appendix [XJ every computed mean value is known exactly and was 
verified to be reproduced by the Monte Carlo simulations within errors. We do 
not plot results for K. They are too boring: errors not visible on the graph and 
exact results agreeing with the data within 1 and occasionally up to around 2 
sigma. 

For the algorithm, the non-interacting case does not seem to be fundamentally 
different from the interacting Gross-Neveu model. All details necessary for this 
extension are given, but the numerical implementation is deferred to a future 
investigation. What remains to be seen is how the correlation between monomers 
of different flavor influences the Monte Carlo dynamics at stronger coupling. 
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5.1 Critical slowing 



We performed a series of simulations of one species of free Majorana fermions at the 
critical value m = 0. In this case, the only infrared scale is given by the system size 
T = L = 8 . . . 128. We simulated the trivial ensembles corresponding to the loop 
class £oo- Results are summarized in table [21 Each run with the local algorithm, 
passing through the lattice in lexicographic order, consists of 10® sweeps of which 
a small fractioiifl is discarded for thermalization. The autocorrelation time r^t.^- 
has been defined and measured as described in |16j . 





exact 


local 


cluster 


L 


K 


K 


Tmt,K 


K 


Tint,K 


8 


0.80738 


0.8080(5) 


5.13(7) 


0.80774(32) 


1.55(2) 


16 


0.78914 


0.7885(4) 


14.4(3) 


0.78932(21) 


1.90(2) 


32 


0.77951 


0.7795(4) 


44.0(1.6) 


0.77949(13) 


2.29(3) 


64 


0.77466 


0.7742(5) 


187(17) 


0.77464(8) 


2.84(4) 


128 


0.77223 


0.7721(4) 


444(58) 


0.77221(5) 


3.50(5) 



Table 2: Monomer density and its integrated autocorrelation time for local and 
cluster simulations at m = and lattice sizes T = L. 

As one expects for local algorithms we see a steeply rising autocorrelation time 
hinting at a dynamical exponent not too far from two — we have no ambition here 
to determine it precisely which would be very costly. One notices, that the error 
(at fixed sweep number) is almost independent of L. This means the variance just 
compensates the growing autocorrelation time and decays roughly proportionally 
to 1/TL. This is in fact implied by scaling and the canonical dimension of the 
(connected) 2-point function of the scalar density. Although the integral over the 
autocorrelation function at L = 128 does not look too unconvincing, one may 
suspect that our number for Ti^t,K may only be a lower bound for this case. 

The cluster simulations in the last columns consist of 0.6 x 10® sweeps which 
in our implementation takes about the same time as the local runs on a single 
PC. We see small slowly rising autocorrelation times. From the two largest lat- 
tices one would estimate an effective dynamical exponent Zcs ~ 0.30, which is a 
typical value for cluster algorithms. In total only about 15 CPU hours went into 
these demonstrations. All codes have been programmed in MATLAB and the update 
routine has about 100 lines (50 without fluctuating boundary conditions). 

The cluster algorithm performs very similarly also in the other topological 
sectors. We did the T = L = 128 runs with spin boundary conditions = (0, 1) 
(£io) and found Ti^t,K = 2.68(4) and = (1, 1) (Cu) with Tint,K = 2.29(3), even 
shorter than the trivial sector. 

■^We discard the first (T/16)^ x 1000 sweeps in the local runs. 
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The next series of runs to be reported is on T = L = 128 lattices at several 
positive and negative masses. Again each data point is produced by 0.6 x 10^ 
sweeps. These simulations included fluctuating boundary conditions. Recorded 
observables were the monomer density K, the boundary conditions and the 
topological flags fap- From these data the distribution of boundary conditions can 
be deduced and we checked their correctness. As an example we consider here 
f l3.2ip and compute for the right hand side 



where on the right hand side the Ising ensemble with fluctuating boundary con- 
ditions e is meant. The result at m = is Siq = 0.76987(6) with the exact value 
being 0.769800361... Errors for this combination of observables are estimated as 
discussed in [16], where the deflnition of Tint,5io from the fluctuations relevant for 
this quantity can be found. We show these autocorrelation times in flgure [3l 
There seems to be a steep rise by about one unit close to m = 0. The second plot 
shows a better resolution of its vicinity. All these numbers stay comfortably small. 
The combination measured here has no serious sign problem. One could however 
construct positive cluster estimators for numerator and denominator. 

A simple example for the use of an improved estimator exploiting the topolog- 
ical information can be given by the two observables with equal mean 



The fractions on the right hand are the part of admissible boundary conditions 
given hj e = (1, 0). Flag positions not tagged here do not allow this value at all. 
At m = we flnd in our run for the left hand side 0.18650(66) [Tint = 0.852(6)] 
and for the other mean 0.18650(37) [T^^t = 1.38(1)] while the exact answer is 
0.186455866.... Of course, the two estimates are strongly correlated. It is clear that 
the left estimate is 'more stochastic' using the actually picked boundary conditions 
in the run and Tint is hence smaller. This pattern is typical for cluster estimators 
with the reduced variance usually overcompensating this effect. 

5.2 Four fermion interaction 

To simulate the interacting Gross-Neveu model (11. ip we represent each factor Z| 
in (11.51) by spin-ensembles. Their e are summed over independently with weights 
depending on the desired boundary conditions for the fermions. At each site x 
an overall conflguration holds ^ ?7,(a;) ^ iV monomers. The only thing that the 




(5.1) 




(5.2) 
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Figure 3: Integrated autocorrelation time for the scalar density (15. ip corresponding 
to fermion boundary conditions antiperiodic in one direction and periodic in the 
other at T = L = 128. 



interaction does is to change the corresponding weight into a factor 



exp 



2 dm'^ 



(2 + m) 



[n/2] 

E 

j=0 



.n-2j 



c{n,m,g). (5.3) 



We update one of the spin flavors at a time. Let us introduce the number n{x) of 
monomers in the other momentarily frozen spins. Then in our update the effective 
monomer weight (or local mass) has to be taken as 



The first few terms are 



(p\o = 2 + m, ip\i = 2 + m + 



9' 



2 + m' 



c{n{x) + 1, m, g) 
c(ji{x),m, g 



ip\2 = 2 + m + 2g 



2 2 + m 
(2 + m) + g^'' 



(5.4) 



(5.5) 



In interacting theories the mass of Wilson fermions undergoes additive renormal- 
ization. Both in the Thirring model {N = 2) and in the higher N Gross-Neveu 
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model perturbative and nonperturbative calculations [1] as well as large ap- 
proximations yield negative values for the critical mass mc{g'^) close to which one 
wants to simulate. This is why it was crucial to extend the cluster algorithm to 
also accommodate ip <2. 

5.3 Additional gauge coupling 

An amusing exercise is to add to the Majorana 'flavor' group a 'color' Z(2) gauge 
interaction. We describe it here only very briefly. We now make use of the gauge 
links in (12. ip and consider 

^coior= ^^^^^'{Zl[m,T]r (5.6) 

with the plaquette fleld Tp made from the links T{x,fi) = ±1. The self-interaction 
can always be added, changing only the monomer weights. We suppress it here. 
Each dimer loop receives as an additional factor the Wilson loop made of r. By 
Stokes' theorem they can be replaced by a product of Tp either on all plaquettes 
where a + spin resides or on those where the negative spins sit, see flgure[Tl Let 
us choose +. In the case of antiperiodic boundary conditions additional loops 
around the torus appear where the torus closes (where the clamps were). As a 
two-dimensional gauge theory is rather trivial the r sum can be carried out. We 
need to know how many plaquettes get tiled with Tp an odd number of times. Let 
us introduce the 'composite' flavor spin- fleld 

S{x) = s(i)(x)s(2)(x) ■ ■ ■ s(^)(a;). (5.7) 

Then 

M± = J2^six),± (5.8) 

X 

contain this information and the extra weight from summing over r is 

22^^[cosh(/?)^+ smh{(3)^- + cosh(/?)*^- smh{(3)^'+]. 
This is equivalent to a magnetic fleld 

if = -^lntanh(/3) (5.9) 

coupled to 5* and fluctuating in sign. In addition we only get contributions if an 
even number of flavors is antiperiodic in each direction separately (due to Z(2) 
conflnement). The fluctuating magnetic fleld can presumably be included in the 
cluster simulation without problems. As we update a given flavor, each spin can 
in addition bond to one exterior ('phantom') spin. 
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6 Conclusion and Outlook 



It seems that with the new cluster algorithm the Gross-Neveu model (and the 
Thirring model) are wide open for high precision simulation. Further observables, 
in particular cluster estimators, remain to be constructed. The O(A^) Noether cur- 
rents may be accessible, for example. Also ratios of partition functions like Z^^/Z}^ 
etc. are expected to have a continuum limit and may serve as renormalized finite 
volume couplings. Majorana fermions are prominent in supersymmetry. Maybe 
some studies on two dimensional supersymmetry become possible with simulations 
very close to the continuum limit. 

An obvious question that every reader will have is if any of this carries over to 
higher dimension and/or more complicated gauge interactions. Let us first caution 
here: fermions in one space dimension are very special. This has manifested itself 
in other so-called Fermi-Bose equivalences in two dimensions. At the heart of 
this in operator language is the fact that the Jordan- Wigner [T7] transformation 
transforms anticommuting to commuting degrees of freedom without generating 
non-localities. This has no obvious generalization to higher dimension (see however 
|18]). The fact that we find positive weights for all topologically trivial loops in 
a way seems to be the euclidean counterpart. Also that Majorana fermions are 
in some sense equivalent to Ising spins it not new, of course, see [12], [SO]- The 
main achievement here is that we simulate in a standard (lattice) euclidean fermion 
formulation and can get really critical. In higher dimension a dimer representation 
for fermions can probably be constructed along similar lines, but the weight will 
be sign-fluctuating in a more essential fashion. The slight hope may be that one 
could be able to handle this sign problem with cluster estimators. In addition, the 
coupling to gauge fields contributing fluctuating Wilson loops is an open problem. 
There are ongoing efforts to simulate discrete models dual to nonabelian gauge 
theories (spin foam) (see ^21j and references therein). This may be an interesting 
view on gauge theories in the context of the approach to fermions developed here. 
In any case, the goal is attractive enough to warrant further thought. 

We would like to acknowledge discussions about the Gross-Neveu model on 
the lattice with Francesco Knechtli, Bjorn Leder, Rainer Sommer and most of all 
Tomasz Korzec. We also thank the Deutsche Forschungsgemeinschaft (DFG) for 
support in the framework of SFB Transregio 9. 

A Exact Majorana partition functions 

The Pfaffian is defined for an antisymmetric matrix Aij of even size, i,j = 1, • • • , 2n, 
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where we integrate over 2n Grassmann variables and a sign convention for the 
measure is imphed. By a change of variables ^ FC, one sees the well known 
identity 

Pf(F^AF) = det(F) Pf(A), (A.2) 
and by squaring the integral in flA.ip . 

[Pf(A)]2 = det(v4) (A.3) 

follows. 

We are interested in A = C(7^(9^ + m — ^d*d) with boundary conditions e. The 
determinant immediately follows by Fourier diagonalization 

detiA) = l[{f + M{pf) (A.4) 



with 

p = — (no + £0/2), — (ni +£1/2) ) , no = 0, . . . ,T - 1, rii = 0, . . . , L - 1 



(A.5) 



and 



Pf, = sin(p^), M{p) =m + ^p^, = 2 sin(p^/2). (A.6) 

In this article we also want to know the relative signs of Pf (A) for different bound- 
ary conditions and possibly negative m. This information is lost in the determinant 
and we proceed differently. 

We define the unitary Fourier transformation matrix 

= ^ e^^^ (A.7) 

which is augmented trivially to include a unit matrix in spin space. Then, using 
( [Ql) we get 

det(F) Pf{A) = Pf (i) (A.8) 

with 

= C(7^p^ + M{p))Sp+g,o- (A.9) 

The factor det(F) is just a phase that we fix later. 

The matrix A consists of antisymmetric blocks where momenta p, —p get paired. 
If they are different (modulo 2tt) such a block contributes one factor (p^ + M(p)^) 
to the Pfaffian. 
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For simplicity we now restrict our discussion to both T and L even. Then all 
momenta get paired non-trivially except p = (0,0), (vr,0), (0,7r), (vr, vr) in the all 
periodic case = (0, 0). We may summarize this result by 

pf.(^) = n vftmw X { f + ■ (a.io) 

Possible extra phase factors can be excluded by considering the limit m — oo. 
Note that the double periodic Pfaffian changes sign at m = 0. 

In this article we need for comparison the scalar condensate in ensembles that 
in fermion language have the partition function 



Z^[m] = Y,c{e)Zl[m]. (A.ll) 



It is given by 



with 



1 / T ^ Id - iy^c{e)z{e,m){CCi) 



z{e,m)= ^'J (A.13) 



and 



p 

All these exact results can be easily evaluated for any finite lattice that is simu- 
lated. For £^ = (0, 0), m — >• the product z remains finite but requires 
precaution numerically. 

The continuum limit of z can presumably be computed analytically. We here 
content ourselves with the numerical construction of their Symanzik expansion in 
a few cases. We set T = L (aspect ratio one) and 

z{e, m) ^ ^ 4(/t, £)L"^ (A.15) 

fc>0 

and compile some values in table El 

The analysis of the asymptotic series was carried out as described in appendix D 
in [22]. The e missing in the table can be computed from those given. Digits are 
quoted such that there is at most an uncertainty of one in the last digit. For k = 
odd corrections vanish. 
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K 


5 













3135575596 


-0 219897 


-1.1 


K 




(in 


di 


do 


1 




12681663 


-0 048092 


-0.108 


1 


(1.0) 


0.27632955 


0.012935 


-0.135 


-1 


(0,0) 


-0.16991195 


-0.08631 


0.183 


-1 


(1,0) 


0.37023294 


0.030368 


-0.285 



Table 3: Symanzik expansion coefficients for ratios of partition functions at differ- 
ent boundary conditions in the finite volume continuum limit at fixed k — mL. 
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